We start by preparing the data from the outputs of scripts 7.1 and 8.1. The results table shown as output is the key to read the figures.
# Results dir + mun shapefile
mun <- st_read("../data/mun/munic_SHP_clean.shp", quiet = TRUE)
# Classes
classes <- read_csv("../config/rcl_tables/land_use/recode_table.csv")
forest_classes <- read_csv("../config/rcl_tables/land_use/recode_table_forest.csv") %>%
rename(new_code = "ID", new_class = "Name")
classes_unique <- unique(classes[, c("new_code", "new_class")]) %>%
bind_rows(forest_classes)
# results key
sce_code_vec <-
as.vector(t(outer(c("BAU-", "R-", "CorPr-", "R-CorPr-", "R(T)-CorPr-"),
c("Hist", "Base", "RCP8"), paste0)))
sce_code_vec_run <-
rep(c("BAU", "R", "CorPr", "R-CorPr", "R(T)-CorPr"), each=3)
results_clean <- read_rds("../data/temp/stsim_run_results.RDS") %>%
dplyr::select(scenarioId, name) %>%
mutate(name = gsub(.$name, pattern = " \\(.+\\)", replacement = "")) %>%
mutate(sce = paste0("sce_", .$scenarioId)) %>%
mutate(chapter = c("none", rep("both", 3), rep("chap_1", 3), rep("chap_2", 9))) %>%
mutate(splitted = str_split(name, " \\| ")) %>%
mutate(climate = unlist(map(splitted, ~unlist(.x[2])))) %>%
mutate(run = unlist(map(splitted, ~unlist(.x[1])))) %>%
dplyr::select(-splitted) %>%
replace_na(list(climate = "none")) %>%
mutate(code = c("Control", sce_code_vec)) %>%
mutate(code_run = c("control", sce_code_vec_run))
# Final extracted datasets
df_final <- readRDS("../outputs/final/final_df_current_density.RDS") %>%
mutate(timestep = (timestep*10)+1990, source = "model")
df_final_origin <- readRDS("../outputs/final/final_df_origin_current_density.RDS") %>%
mutate(timestep = timestep*10+1990, source = "model")
# Stratum key
stratum_key <- read_csv("../config/stsim/SecondaryStratum.csv") %>%
mutate(ID = as.factor(ID)) %>%
rename(zone=ID, MUS_NM_MUN=Name) %>%
left_join(df_final, by="zone") %>%
filter(MUS_NM_MUN!="Not Monteregie")
# Sum all municipalities
df_summarised <- df_final %>%
group_by(sce, timestep, species, iteration) %>%
summarise(sum_cur = sum(current)) %>% ungroup %>%
mutate(source = "model")
df_origin_summarised <- df_final_origin %>%
group_by(timestep, species) %>%
summarise(sum_cur = sum(mean)) %>% ungroup %>%
mutate(sce = "sce_0", source = "observation")
# Joined dataset
joined <- full_join(df_summarised, df_origin_summarised,
by=c("sce", "source", "species", "timestep", "sum_cur")) %>%
left_join(results_clean, by = "sce") %>%
replace_na(list(climate = "none", run = "historic run")) %>%
# Make factors
mutate(sce = as.factor(sce), run = as.factor(run), sum_cur = 10*(sum_cur),
climate = factor(climate, levels = c("none", "historic",
"baseline", "RCP 8.5")),
run = factor(run, levels = c("historic run", "BAU run",
"BAU run + corrs protection", "BAU run + ref",
"BAU run + corrs protection + ref",
"BAU run + corrs protection + ref TARGETED")))
# Roc curves data
urb <- readRDS("../data/temp/fit_rs_outcome_rf_urb_2.RDS")
agex <- readRDS("../data/temp/fit_rs_outcome_rf_agex_2.RDS")
# Histogram data
hist_original <- read_csv("../outputs/final/final_values_output_original.csv")
hist_true <- read_csv("../outputs/final/final_values_output_TRUE.csv") %>%
mutate(sce = "TRUE")
histograms <- read_csv("../outputs/final/final_values_output.csv") %>%
bind_rows(hist_original) %>%
bind_rows(hist_true) %>%
left_join(results_clean, by="sce")
# SURF Data
surf <- read_csv("../surf/surf_output.csv") %>%
mutate(timestep = timestep*10+1990) %>%
left_join(results_clean, by=c("scenario"="scenarioId")) %>%
replace_na(list(climate = "none", run = "historic run")) %>%
# Make factors
mutate(sce = as.factor(sce), run = as.factor(run),
climate = factor(climate, levels = c("none", "historic",
"baseline", "RCP 8.5")),
run = factor(run, levels = c("historic run", "BAU run",
"BAU run + corrs protection", "BAU run + ref",
"BAU run + corrs protection + ref",
"BAU run + corrs protection + ref TARGETED")))
# Bar plot data
bar_data <- read_csv("../outputs/final/final_bar_plot_data.csv")
results_clean
the_width = 18
the_height = 15
the_dpi = 500
fig_1_historic <- joined %>%
filter(climate == "none") %>%
group_by(scenarioId, species, iteration) %>%
group_modify(~mutate(.x , baseline = subset(.x, timestep == 1990)$sum_cur,
the_diff = ((sum_cur-baseline)/baseline)*100 )) %>%
ggplot(aes(x=timestep, y=the_diff, col=source)) +
scale_color_manual(values=c('#d8b365','#5ab4ac'),
labels = c("Model", "Observation")) +
geom_smooth(aes(group = sce), method = "glm", se=FALSE) +
scale_linetype_manual(values = c(1,3,2,5))+
geom_point(aes(group = sce)) +
facet_wrap(~species, scales = "fixed") +
labs(y = "Current flow (% Change)",
x = "Year",
col = "Source") +
NULL
fig_1_historic
Connectivity change for species through time, 1990-2010
ggsave(fig_1_historic,
filename = "../thesis/figures/connectivity_decrease_x5species_historic.png",
width = the_width, height = the_height, dpi = the_dpi)
colors_sce <-
data.frame(sce = c("Hist", "Baseline", "RCP85"),
color = c("#8DA0CB", "#FC4E2A", "#800026"), stringsAsFactors = F)
fig_1_static_chap_1 <- joined %>%
filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>%
group_by(scenarioId, species, iteration) %>%
group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur,
the_diff = ((sum_cur-baseline)/baseline)*100 )) %>%
ggplot(aes(x=timestep, y=the_diff, col=climate)) +
scale_color_manual(values = colors_sce$color,
labels = c("Historic", "Baseline", "RCP 8.5")) +
geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
#geom_point(aes(group = sce)) +
scale_linetype_manual(values=c("solid", "twodash")) +
facet_wrap(~species, scales = "fixed") +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
#geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
labs(y = "Current flow (% Change)",
x = "Year",
col = "Climate",
linetype = "Run") +
guides(col = guide_legend(ncol = 2, nrow=2),
linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
theme(legend.position = c(0.85, 0.15))+
NULL
fig_1_static_chap_1
Connectivity change for species through time, 2010-2100
ggsave(fig_1_static_chap_1,
filename = "../thesis/figures/connectivity_decrease_x5species_chap1.png",
width = the_width, height = the_height, dpi = the_dpi)
# Default size
gridline.label.offset = 4
legend.text.size = 20
group.line.width = 0.9
grid.label.size = 10
background.circle.colour = "#ffffff"
group.point.size = 5
axis.label.size = 8
group.colours = RColorBrewer::brewer.pal(name = "Paired",n=10)[c(2,4,6,8,10)]
radar_data_chap1 <- joined %>%
filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>%
group_by(timestep, species, code) %>%
summarise(sum_cur = mean(sum_cur)) %>% ungroup %>%
pivot_wider(names_from = timestep, values_from = sum_cur) %>%
mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>%
select(-c(`2010`:`2100`)) %>%
pivot_wider(names_from = code, values_from = diff) %>%
rename(group = species)
radar_chap1 <-
ggradar(radar_data_chap1, centre.y = -20, legend.position = "right",
grid.min = -20, grid.max = 3, grid.mid = 0,
values.radar = c("-20 %", "0 %", ""),
group.colours = group.colours,
gridline.label.offset = gridline.label.offset,
legend.text.size = legend.text.size,
group.line.width = group.line.width,
grid.label.size = grid.label.size,
background.circle.colour = background.circle.colour,
group.point.size = group.point.size,
axis.label.size = axis.label.size)
radar_chap1
ggsave("../thesis/figures/radar_ggradar_chap1.png", radar_chap1,
width = the_width, height = the_height, dpi = the_dpi)
fig_1_static_chap_2 <- joined %>%
filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>%
group_by(scenarioId, species, iteration) %>%
group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur,
the_diff = ((sum_cur-baseline)/baseline)*100 )) %>%
ggplot(aes(x=timestep, y=the_diff, col=climate)) +
scale_color_manual(values = colors_sce$color,
labels = c("Historic", "Baseline", "RCP 8.5")) +
geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
facet_wrap(~species, scales = "fixed") +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
labs(y = "Current flow (% Change)",
x = "Year",
col = "Climate",
linetype = "Run") +
guides(col = guide_legend(ncol = 2, nrow=2),
linetype = guide_legend(nrow = 2, ncol=2,
override.aes = list(colour = 'black'))) +
theme(legend.position = c(0.85, 0.15))+
NULL
fig_1_static_chap_2
Connectivity change for species through time, 2010-2100
ggsave(fig_1_static_chap_2,
filename = "../thesis/figures/connectivity_decrease_x5species_chap2.png",
width = the_width, height = the_height, dpi = the_dpi)
radar_data_chap2 <- joined %>%
filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>%
group_by(timestep, species, code) %>%
summarise(sum_cur = mean(sum_cur)) %>% ungroup %>%
pivot_wider(names_from = timestep, values_from = sum_cur) %>%
mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>%
select(-c(`2010`:`2100`)) %>%
pivot_wider(names_from = code, values_from = diff) %>%
rename(group = species)
radar_chap2 <-
ggradar(radar_data_chap2, centre.y = -20, legend.position = "right",
grid.min = -20, grid.max = 3, grid.mid = 0,
values.radar = c("-20 %", "0 %", ""),
group.colours = group.colours,
gridline.label.offset = gridline.label.offset,
legend.text.size = legend.text.size,
group.line.width = group.line.width,
grid.label.size = grid.label.size,
background.circle.colour = background.circle.colour,
group.point.size = group.point.size,
axis.label.size = axis.label.size-2)
radar_chap2
ggsave("../thesis/figures/radar_ggradar_chap2.png", radar_chap2,
width = the_width, height = the_height, dpi = the_dpi)
radar_data_all <- joined %>%
filter(climate != "none") %>%
group_by(timestep, species, code) %>%
summarise(sum_cur = mean(sum_cur)) %>% ungroup %>%
pivot_wider(names_from = timestep, values_from = sum_cur) %>%
mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>%
select(-c(`2010`:`2100`)) %>%
pivot_wider(names_from = code, values_from = diff) %>%
rename(group = species)
radar_all <-
ggradar(radar_data_all, centre.y = -20, legend.position = "right",
grid.min = -20, grid.max = 3, grid.mid = 0,
values.radar = c("-20 %", "0 %", ""),
group.colours = group.colours,
gridline.label.offset = gridline.label.offset,
legend.text.size = legend.text.size,
group.line.width = group.line.width,
grid.label.size = grid.label.size,
background.circle.colour = background.circle.colour,
group.point.size = group.point.size,
axis.label.size = axis.label.size-2)
radar_all
ggsave("../thesis/figures/radar_ggradar_both.png", radar_all,
width = the_width, height = the_height, dpi = the_dpi)
ROC curves from fitting both models, plotting with patchwork package.
p1 <- bind_rows(urb, .id = "fold") %>%
mutate(Fold = factor(fold, levels = as.character(1:10))) %>%
group_by(Fold) %>%
roc_curve(.pred, truth = outcome_fact) %>%
autoplot(add=T) +
ggtitle(paste("Urbanisation")) +
annotate(x = 0.75, y = 0.25, geom="label", size = 3,
label = as.character(paste("Av AUC = 0.938 +/- 0.002"))) +
theme(legend.position = "none") +
NULL
p2 <- bind_rows(agex, .id = "fold") %>%
mutate(Fold = factor(fold, levels = as.character(1:10))) %>%
group_by(Fold) %>%
roc_curve(.pred, truth = outcome_fact) %>%
autoplot(add=T) +
ggtitle(paste("Agricultural Expansion")) +
annotate(x = 0.75, y = 0.25, geom="label", size = 3,
label = as.character(paste("Av AUC = 0.929 +/- 0.002"))) +
NULL
full = p1 + p2
full
ggsave("../thesis/figures/double_roc_resample.png", full)
bi_colors <- c("#67a9cf", "#ef8a62")
hist_plot_origin <- histograms %>%
filter(sce %in% c("TRUE", "sce_37")) %>%
filter(ts == 2) %>%
mutate(ts = ts*10+1990) %>%
mutate(ts = as.factor(ts)) %>%
ggplot(aes(x=bins, y=n, fill=sce, colour =sce)) +
geom_density(stat='identity', show.legend=T, alpha=0.3)+
scale_fill_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
scale_color_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
facet_wrap(~species) +
labs(x = "Flow intensity distribution (log)",
y = "") +
theme(strip.text.y = element_text(angle=360, size=10, hjust = 0),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y.left = element_blank(),
axis.ticks.y.left = element_blank()) +
NULL
hist_plot_origin
ggsave("../thesis/figures/original_hists.png", hist_plot_origin,
width = the_width, height = the_height, dpi = the_dpi)
bi_colors <- c("#67a9cf", "#ef8a62")
hist_plot_1 <- histograms %>%
mutate(ts = as.factor(ts)) %>%
filter(chapter %in% c('both', 'chap_1')) %>%
ggplot(aes(x=bins, y=n, group=ts)) +
geom_density(stat='identity', show.legend=T, aes(fill=factor(ts), color=factor(ts)), alpha=0.3)+
facet_grid(code~species) +
scale_fill_manual(values = bi_colors,
labels = c("2010", "2100")) +
scale_color_manual(values = bi_colors,
labels = c("2010", "2100")) +
labs(x = "Flow intensity distribution (log)",
y = "")+
theme(strip.text.y = element_text(angle=360, size=10, hjust = 0),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y.left = element_blank(),
axis.ticks.y.left = element_blank(),
legend.position = c(0.95, 0.96)) +
labs(color="Timestep", fill="Timestep") +
NULL
hist_plot_1
ggsave("../thesis/figures/hist_chap1.png", hist_plot_1, width = 18, height = 15)
hist_plot_2 <- histograms %>%
mutate(ts = as.factor(ts)) %>%
filter(chapter %in% c('chap_2')) %>%
ggplot(aes(x=bins, y=n, group=ts)) +
geom_density(stat='identity', aes(fill=factor(ts), colour=factor(ts)), alpha=0.3)+
facet_grid(code~species) +
scale_fill_manual(values = bi_colors,
labels = c("2010", "2100")) +
scale_color_manual(values = bi_colors,
labels = c("2010", "2100")) +
labs(x = "Flow intensity distribution (log)",
y = "") +
theme(strip.text.y = element_text(angle=360, size=10, hjust = 0),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y.left = element_blank(),
axis.ticks.y.left = element_blank(),
legend.position = c(0.95, 0.96)) +
labs(color="Timestep", fill="Timestep") +
NULL
hist_plot_2
ggsave("../thesis/figures/hist_chap2.png", hist_plot_2, width = 18, height = 20)
surf_1 <- surf %>%
mutate(scenario = as.factor(scenario)) %>%
filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>%
filter(climate != "none") %>%
group_by(scenario, species, iter) %>%
group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb,
the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>%
ggplot(aes(x = timestep, y = the_diff, color=climate)) +
scale_linetype_manual(values=c("solid", "twodash")) +
geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
#geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
facet_wrap(~species, scales = "fixed")+
scale_color_manual(values = colors_sce$color,
labels = c("Historic", "Baseline", "RCP 8.5")) +
labs(y = "Change in nb of Keypoints detected (%)",
x = "Year",
col = "Climate",
linetype = "Run") +
guides(col = guide_legend(ncol = 2, nrow=2),
linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
NULL
surf_1
ggsave("../thesis/figures/surf_chap1.png", surf_1,
width = the_width, height = the_height, dpi = the_dpi)
surf_radar_data_1 <- surf %>%
filter(climate != "none", chapter %in% c("none", "both","chap_1")) %>%
filter(sce != "sce_0", name != "historic run") %>%
group_by(timestep, species, code) %>%
summarise(kp_nb = mean(kp_nb)) %>% ungroup %>%
pivot_wider(names_from = timestep, values_from = kp_nb) %>%
mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>%
select(-c(`2010`:`2100`)) %>%
pivot_wider(names_from = code, values_from = diff) %>%
rename(group = species)
surf_radar_1 <-
ggradar(surf_radar_data_1, centre.y = -60, legend.position = "right",
grid.min = -60, grid.max = 5, grid.mid = 0,
values.radar = c("-60 %", "0 %", ""),
group.colours = group.colours,
gridline.label.offset = 20,
legend.text.size = legend.text.size,
group.line.width = group.line.width,
grid.label.size = grid.label.size,
background.circle.colour = background.circle.colour,
group.point.size = group.point.size,
axis.label.size = axis.label.size)
surf_radar_1
ggsave("../thesis/figures/surf_radar_chap1.png", surf_radar_1,
width = the_width, height = the_height, dpi = the_dpi)
surf_2 <- surf %>%
mutate(scenario = as.factor(scenario)) %>%
filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>%
filter(climate != "none") %>%
group_by(scenario, species, iter) %>%
group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb,
the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>%
ggplot(aes(x = timestep, y = the_diff, color=climate)) +
geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
#geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
facet_wrap(~species, scales = "fixed")+
scale_color_manual(values = colors_sce$color,
labels = c("Historic", "Baseline", "RCP 8.5")) +
labs(y = "Change in nb of Keypoints detected (%)",
x = "Year",
col = "Climate",
linetype = "Run") +
guides(col = guide_legend(ncol = 2, nrow=2),
linetype = guide_legend(nrow=2, ncol=2,
override.aes = list(colour = 'black'))) +
NULL
surf_2
ggsave("../thesis/figures/surf_chap2.png", surf_2,
width = the_width, height = the_height, dpi = the_dpi)
surf_radar_data_2 <- surf %>%
filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>%
filter(sce != "sce_0", name != "historic run") %>%
group_by(timestep, species, code) %>%
summarise(kp_nb = mean(kp_nb)) %>% ungroup %>%
pivot_wider(names_from = timestep, values_from = kp_nb) %>%
mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>%
select(-c(`2010`:`2100`)) %>%
pivot_wider(names_from = code, values_from = diff) %>%
rename(group = species)
surf_radar_2 <-
ggradar(surf_radar_data_2, centre.y = -60, legend.position = "right",
grid.min = -60, grid.max = 5, grid.mid = 0,
values.radar = c("-60 %", "0 %", ""),
group.colours = group.colours,
gridline.label.offset = 22,
legend.text.size = legend.text.size,
group.line.width = group.line.width,
grid.label.size = grid.label.size,
background.circle.colour = background.circle.colour,
group.point.size = group.point.size,
axis.label.size = axis.label.size-2)
surf_radar_2
ggsave("../thesis/figures/surf_radar_chap2.png", surf_radar_2,
width = the_width, height = the_height, dpi = the_dpi)
source('~/Documents/Master/Thesis/land_con_monteregie/scripts/functions/draw_scenario.R')
for (scenario in c(39, 42, 45, 48, 51)){
the_plot <- draw_scenario(scenario, mxp = mxp)
the_name <- paste0("scenario_", as.character(scenario),"_compare.png")
ggsave(plot = the_plot,
filename = paste0("../thesis/figures/",the_name), width = 20, height = 10, dpi = the_dpi)
}
## Parsed with column specification:
## cols(
## source = col_character(),
## name = col_character(),
## class = col_character(),
## subclass = col_logical(),
## old_code = col_double(),
## new_code = col_double(),
## new_class = col_character()
## )
## Parsed with column specification:
## cols(
## ID = col_double(),
## Name = col_character()
## )
## Parsed with column specification:
## cols(
## source = col_character(),
## name = col_character(),
## class = col_character(),
## subclass = col_logical(),
## old_code = col_double(),
## new_code = col_double(),
## new_class = col_character()
## )
## Parsed with column specification:
## cols(
## ID = col_double(),
## Name = col_character()
## )
## Parsed with column specification:
## cols(
## source = col_character(),
## name = col_character(),
## class = col_character(),
## subclass = col_logical(),
## old_code = col_double(),
## new_code = col_double(),
## new_class = col_character()
## )
## Parsed with column specification:
## cols(
## ID = col_double(),
## Name = col_character()
## )
## Parsed with column specification:
## cols(
## source = col_character(),
## name = col_character(),
## class = col_character(),
## subclass = col_logical(),
## old_code = col_double(),
## new_code = col_double(),
## new_class = col_character()
## )
## Parsed with column specification:
## cols(
## ID = col_double(),
## Name = col_character()
## )
## Parsed with column specification:
## cols(
## source = col_character(),
## name = col_character(),
## class = col_character(),
## subclass = col_logical(),
## old_code = col_double(),
## new_code = col_double(),
## new_class = col_character()
## )
## Parsed with column specification:
## cols(
## ID = col_double(),
## Name = col_character()
## )
final_df_39 <- bar_data %>%
filter(scenario == 39)
bar_plot_39 <- final_df_39 %>%
#filter(value > 10) %>%
mutate(value = as.factor(value), timestep = as.factor(timestep)) %>%
ggplot(aes(fill=timestep, y=area, x=new_class)) +
geom_bar(position="dodge", stat="identity") +
#geom_errorbar(aes(ymin=mean_count-sd_count, ymax=mean_count+sd_count), width=.2,
# position=position_dodge(.9)) +
scale_fill_manual(values = bi_colors) +
theme(legend.position = c(0.85, 0.7)) +
labs(x = "Land use Class",
y = "Mean Area (hectare)",
fill = "Year") +
NULL
bar_plot_39
#ggsave("../thesis/figures/bar_BAU.png", bar_plot_39, width = 15, height = 10)
# draw_scenario(42)
final_df_42 <- bar_data %>%
filter(scenario == 42)
bar_plot_42 <- final_df_42 %>%
#filter(value > 10) %>%
mutate(value = as.factor(value), timestep = as.factor(timestep)) %>%
ggplot(aes(fill=timestep, y=area, x=new_class)) +
geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values = bi_colors) +
#geom_errorbar(aes(ymin=mean_count-sd_count, ymax=mean_count+sd_count), width=.2,
# position=position_dodge(.9)) +
theme(legend.position = c(0.85, 0.7)) +
labs(x = "Land use Class",
y = "Mean Area (hectare)",
fill = "Year") +
NULL
bar_plot_42
#ggsave("../thesis/figures/bar_Ref.png", bar_plot_42, width = 15, height = 10)
final_df_39 <- final_df_39 %>%
mutate(sce = "BAU") %>%
mutate(type = ifelse(value <10, "non_forest", "forest"))
final_df_42 <- final_df_42 %>%
mutate(sce = "Reforestation") %>%
mutate(type = ifelse(value <10, "non_forest", "forest"))
final_all <- bind_rows(final_df_39, final_df_42)
forest <- final_all %>%
filter(type == "forest") %>%
mutate(value = as.factor(value), timestep = as.factor(timestep)) %>%
ggplot(aes(fill=timestep, y=area, x=new_class)) +
geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values = bi_colors) +
theme(legend.position = "none") +
facet_grid(~sce) +
labs(x = "Land use Class",
y = "Mean Area (hectare)",
fill = "Year") +
NULL
non_forest <- final_all %>%
filter(type != "forest") %>%
mutate(value = as.factor(value), timestep = as.factor(timestep)) %>%
ggplot(aes(fill=timestep, y=area, x=new_class)) +
geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values = bi_colors) +
theme(legend.position = "right") +
facet_grid(~sce) +
labs(x = "Land use Class",
y = "Mean Area (hectare)",
fill = "Year") +
NULL
assembled <- wrap_plots(forest, non_forest, ncol=1, nrow=2)
assembled
ggsave("../thesis/figures/bar_both.png", assembled,
width = the_width, height = 20, dpi = the_dpi)